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Abstract. In this article we present a parallel modular algorithm to com- 
pute all solutions of a given zero-dimensional polynomial system of equations 
over the rationals. In fact, we compute a triangular decomposition of the 
corresponding ideal in the polynomial ring over the rationals using modular 
methods, and then apply a solver for univariate polynomials. 

1. Introduction 

One possible approach^ to find the solutions of a zero-dimensional system of 
multivariate polynomials is the triangular decomposition of the corresponding ideal 
since triangular systems of polynomials can be solved using a univariate solver 
recursively. Therefore, we recall the definition and some properties of a triangular 
decomposition. For details we refer to |GP07| . |M93| and |M97j . 

Let K he a field, X — {xi,...,Xn} a set of variables, / C K[X] a zero- 
dimensional ideal, and we fix > to be the lexicographical ordering induced by 
Xi > . . . > Xn- If / S K[X] is a polynomial, then we denote by LC(/) the leading 
coefficient of /, by LM(/) the leading monomial of /, andbyLT(/) = LC(/)-LM(/) 
the leading term of /. 

Definition 1.1. A set of polynomials F = {/i, . . . , /„} C K[X] is called a trian- 
gular set if LT(/i) = x^Lj^]^ for some ai > and each i — 1, . . . , n. 

A list T ~ {Fi, . . . , Fs) of triangular sets is called a triangular decomposition of 
a zero-dimensional ideal I C K[X] if Vl = \/ {Fi) n . . . n \/{Fs). 

Remark 1.2. 

(1) A triangular set is a Grobner basis. 

(2) A minimal Grobner basis of a maximal ideal is a triangular set, and the 
primary decomposition of vT is a triangular decomposition of /. 

The following two lemmata are the basis for the algorithm by Moller (cf. [M93] ) 
which avoids primary decomposition. 

Lemma 1.3. Let G = {gi, . . . ,gm} be a reduced Grobner basis of the zero-dimen- 
sional ideal I C if [A] such that LM(gi) < ... < LM(g„). Moreover, let gi — 
J2%o9't3x{ with g[j G K[x2, . . . ,Xn], g'ia, ^ 0, and G' = ■ • ■ >5m-ia„-i}- 

Then G' is a Grobner basis of {gi, . . . ,gm-i) ■ g-m, and we have {G' ,gm) = (C, G). 
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^In Singular (cf. |DGPS12] ) this approach is implemented in the library solve. lib on the 
basis of an univariate Laguerre solver (cf. [RR78. §8.9-8.13]). 
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Proof. The proof can be found in [M93[ Lemma 7] . □ 

Lemma 1.4. Let I C K[X] be a zero- dimensional ideal, and h G -f'^i-'^]- Then the 
following hold. 

(1) VT= vW^n/TTTI. 

(2) dim^ {K[X]/I) = dimx {K[X]/{I, h)) + dini,^ {K[X]I{I : h)). 
Proof. 

(1) The proof is an exercise in jGP071 Exercise 4.5.3]. 

(2) We consider two exact sequences. The first one 

^ (/ : /i)// ^ K[X]/I K[X]/I K[X]/{I, h) 

yields dimA' ((/ : h)/l) = dim^ (^K[X]/{I, h)), and the second one 

— >{I : h)/I — > K[X]/I — > K[X]/{I : h) — ^ 

yields dim/f {K[X]/I) = dim^: ((/ : h)/l) + Au^k {K[X]/{I : h)). Sum- 
marized, we obtain 

dim^ {K[X]/I) = dim^ {K[X]/{I, h)) + dim^ {K[X]/{I : h)). 

□ 

Consequently, for hi, . . . ,hi G we are able to apply Lemma 11.41 inductively 

and obtain 

Vl=^{I,hi) n ^/I:hi 

= v/(/,/ii,/i2> n ^{I,hi) : h2n ■ hi 



= v/(/, hi,..., hi) n l^pl ^/{I,hi,...,h,^i) : h,j . 

Together with Lemma 11.31 we conclude the following corollary. 

Corollary 1.5. With the assumption of Lemma \1.3[ let G' \ G — {hi, . . . ,hi}. 
Then the following hold. 

(1) IfG'\G^(!), then I C (G,/ii) and I C (G) -.hi. 

(2) VT = V(G',.9,n) n (n-=i hi,..., : h,) . 

(3) dim^ {K[X]/I) = J2\^^ dim^ {K[X]/{{G, hi,..., K_i) : K)) . 

With regard to Corollary [T31J2), especially (G, hi, . . . ,hi) — (G', g™) with G' C 
K[x2, . ■ . ,Xn] is predestined for induction since a/ {G',gm) = V {Fi-9m) n . . . n 
\/{F[7gm) if J'' = {Fi, ■ ■ ■ ,F^) is a triangular decomposition of G'. Referring to 
Corollary II. 5f 3). the triangular decomposition obtained by iterating the approach 
of Corollary 11.51 respects the multiplicities of the zeros of /. Therefore the zero- 
sets of different triangular sets are in general not disjoint as the following example 
shows. 

Example 1.6. Let G = {xY',xiX2 +X2,xl^} C Q[xi,X2]. Then we obtain G' = 
{x2°,X2}, G' \ G = {X2}, and the triangular decomposition J" = (^1,^2) of 
(G) with Fi = (G) : xl = {xl,xi + x^) and F2 = (G,xi) = {xlx\^). More- 
over, it holds dimQ (Q[xi,X2]/(G)) = 40, dimQ (Q[a;i, X2]/((G) : x^)) = 7, and 
dimQ {Q[xi,X2]/{G,xl,)) =33. 
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Algorithm [T] shows the algorithm by Moller to compute the triangular decompo- 
sition of a zero-dimensional idealH 



Algorithm 1 Triangular Decomposition (triangM) 
Input: / C if [AT], a zero-dimensional ideal . 

Output: T ~ (Fi. . . . , Fs), a triangular decomposition of / C K[X] such that 

dim^ (if[^]//) = E:=i dim^ iK[X]/{F,)). 
1: compute G = {gi, . . . ,gm}, a reduced Grobner basis of / with respect to the 

lexicographical ordering > such that LM(gi) < . . . < LM(g„i); 
2: compute G' — {g'l, . . . , g'm-i} Q K[x2, ■ ■ ■ , Xn] where gl is the leading coefficient 

of ,gi in {K[x2, . . . ,x„])[a;i]; 
3: F' = triangM ((G')); 
4: ^ = {F' U {gm} I F' e F'}; 
5: for 1 < i < m — 1 do 
6: li g[^ G then 
7: J" J" U triangM ((G) : 
8: G = GU{5a; 
9: return F] 



Note that Algorithm [T] is only based on Grobner basis computations, and does 
not use random elements. Hence, the result is uniquely determined which allows 
modular computations. In the following we fix Algorithm [1] to compute a triangular 
decomposition. 

Remark 1.7. Replacing line 7 in Algorihtm [T] by F = FU triangM ((G) : we 
obtain a disjoint triangular decomposition (i.e. Fi,Fj G F with Fi ^ Fj implies 
(Fi) + (Fj) = K[X]). This decomposition does in general not respect the multiplic- 
ities (i.e. dim^ {K[X]/I) ^ E'=idimK {K[X]/ {F,))). 

2. Modular methods 

One possible modular approach for solving a zero-dimensional ideal is to just 
replace each involved Grobner basis computation by its corresponding modular al- 
gorithm as described in [IPS 11) . Particularly, we can replace line 1 in Algorithm 
[T]by G = modStd(/). In this case it is possible to apply the probabilistic variant 
since we can easily verify in the end by a simple substitution whether the solutions 
obtained from the triangular sets are really solutions of the original ideal. Nev- 
ertheless, we propose to compute the whole triangular decomposition via modular 
methods actually. The verification is then similar by just substituting the obtained 
result into the input polynomials. 

We consider the polynomial ring Q[Ar], fix a global monomial ordering > on 
Q[Ar], and use the following notation: If S" C Q[X] is a set of polynomials, then 
LM(S') := {LM(/) | / G 5} denotes the set of leading monomials of 5'. If / e Q[X] 
is a polynomial, / = (/i, . . . , fr) Q is an ideal, and p is a prime number which 

does not divide any denominator of the coefficients of f,fi,...,fr, then we write 
fp := (/ mod p) e ¥p[X] and Ip ((/i)p, . . . , (/,)p) C ¥p[X]. 



'The corresponding procedure is implemented in Singular in the library triang.lib. 
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In the following, / = (/i, . . . , /r) C Q[X] will be a zero-dimensional ideal. The 
triangular decomposition algorithm (Algorithm [1} applied to / returns a list of 
triangular sets J" = (Fi, . . . , F^) such that a/T = Cl^^i and dimQ (Q[X]//) = 

ELidiniQ(Q[X]/(F,)). 

With respect to modularization, the following lemma obviously holds: 

Lemma 2.1. With notation as above, let p he a sujficiently general prime number. 
Then Ip is zero- dimensional, {{Fi)p, . . . , {Fs)p) is a list of triangular sets, and 

Relying on Lemma |2. 11 the basic idea of the modular triangular decomposition 
is as follows. First, choose a set V of prime numbers at random. Second, compute 
triangular decompositions Fp of Ip for p E V. Third, lift the modular triangular 
sets to triangular sets over Q[X]. 

The lifting process consists of two steps. First, the set TV :— {J-p \ p G V} is 
lifted to with {Fi)^ C {Z/N'Z)[X] and N :— Yipt^-pP applying the Chinese 
remainder algorithm to the coefficients of the polynomials occurring in TV. Second, 
we obtain T with Fi C Q[X] by lifting the modular coefficients occurring in Tp 
to rational coefficients via the Farey rational mapQ. This map is guaranteed to be 
bijective provided that -y/ N/2 is larger than the moduli of all coefficients of elements 
in T with F, Q Q[X]. 

We now define a property of the set of primes V which guarantees that the lifting 
process is feasible and correct. This property is essential for the algorithm. 

Definition 2.2. Let T — {Fi, . . . , Fg) be the triangular decomposition of the ideal 
/ computed by Algorithm [1] 

(1) A prime number p is called lucky for / and T if {{Fi)p, . . . , {Fs)p) is a 
triangular decomposition of Ip. Otherwise p is called unlucky for / and T. 

(2) A set V of lucky primes for / and T is called sujficiently large for / and J- 
if 

JJ^ p > max{2 • |c|^ | c coefficient occuring in J-}. 

From a theoretical point of view, the idea of the algorithm is now as follows: 
Consider a sufficiently large set V of lucky primes for / and T, compute the trian- 
gular decomposition of the Ip, p eV, via Algorithm [TJ and lift the results to the 
triangular decomposition of / as aforementioned. 

From a practical point of view, we face the problem that we do not know in 
advance whether a prime number p is lucky for / and J- . 

To handle this problem, we fix a natural number t and an arbitrary set of primes 
V of cardinality t. Having computed FV, we use the following test to modify V 
such that all primes in V are lucky with high probability: 

deleteUnluckyPrimesTriang: We define an equivalence relation on {FV ,V) 
by (Fp,p) ^ (J-,, q) :^ = #7^ and {LM(Fp) | Fp G Fp} = {LM(F,) | F, G 

FqY). Then the equivalence class of largest cardinality is stored in {FV,V), the 
others are deleted. 



Farey fractions refer to rational reconstruction. A definition of Farey fractions, the Farey 
rational map, and remarks on the required bound on the coefiicients can be found in |KG83| . 
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Since we do not know a priori whether the equivalence class chosen is indeed 
lucky and whether it is sufRciently large for / and T , we proceed in the following 
way. We lift the set TV to T over as described earlier, and test the result 

with another randomly chosen prime number: 

pTestTriang: We randomly choose a prime number p ^ V such that p does not 
divide the numerator and denominator of any coefficient occuring in a polynomial 
in {fi, . . . , fr} or J-. The test returns true if {Fp \ F Cz T) equals the triangu- 
lar decomposition J-p computed by the fixed Algorithmic applied on Ip, and false 
otherwise. 

If pTestTriang returns false, then V is not sufficiently large for / and T or 
the equivalence class of prime numbers chosen was unlucky. In this case, we enlarge 
the set V by t new primes and repeat the whole process. On the other hand, if 
pTestTriang returns true, then we have a triangular decomposition F of I with 
high probability. In this case, we compute the solutions S C C" with multiplici- 
ties of the triangular sets F — {/i, . . . , /„} G as follows. Solve the univariate 
polynomial fi{xi) via a univariate solver counting multiplicities, substitute xi in 
f2{xi,X2) by these solutions of fi{xi), solve the corresponding univariate polyno- 
mial, and continue inductively this way (call this step SOLVeTriang). Finally, we 
verify the result S* C C" by testing whether the original ideal / C Q[X] vanishes 
on S" C C", and whether the sum of the multiplicities equals the Q-dimension of 
Q[X]// (caU this step testZero). 

We summarize modular solving in Algorithm [Sj^ 



Algorithm 2 Modular Solving (modSolve) 
Input: / C Q[Ar], a zero-dimensional ideal. 

Output: S C C", a set of points in C" such that f[P) = for aU / e /, P e S. 

1: choose V, a list of random primes; 
2: FV = 0; 

3: loop 

4: for p e P do 

5: J-p = triaiigM(/p), the triangular decomposition of Ip via Algorithm [1] 

6: FV = TV \J {Fp}\ 

7: {FV,V) = deleteUnluckyPrimesTriang(J"7',P); 

8: lift {TV,V) to T over Q[X] by applying Chinese remainder algorithm and 
Farey rational map; 

9: if pTestTriang(/, J", V) then 

10: S = solveTriang(J"); 

11: if testZero(/, S) then 

12: return S; 

13: enlarge V; 



Remark 2.3. In Algorithm [21 the triangular sets Tp can be computed in parallel. 
Furthermore, we can parallelize the final verification whether / C Q[X] vanishes on 
S c C". 



The corresponding procedures are implemented in Singular in the library modsolve . lib. 
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3. Examples and timings 

In this section we provide examples on which we time the algorithm modSolve 
(cf. Algorithm and its parallel version as opposed to the algorithm solve (the 
procedure solve is implemented in Singular in the library solve. lib and com- 
putes all roots of a zero-dimensional input ideal using triangular sets). Timings 
are conducted by using Singular 3-1-4 on an AMD Opteron 6174 machine with 
48 CPUs, 2.2 GHz, and 128 GB of RAM running the Gentoo Linux operating sys- 
tem. All examples are chosen from The SymbolicData Project (cf. tG12j) and the 
PHCpack test suite (cf. [V99) ). 

Remark 3.1. The parallelization of the modular algorithm is attained via multiple 
processes organized by Singular library code. Consequently, a future aim is to 
enable parallelization in the kernel via multiple threads. 

We choose the following examples to emphasize the superiority of modular solv- 
ing and especially its parallelization: 

Example 3.2. Cyclic_7 (cf. [GT2] , [V99]). 

Example 3.3. Verschelde_aoon5 (cf. jG12] . |V99| ). 

Example 3.4. Verschelde_noon6 (cf. [GT2] . jV99] ). 

Example 3.5. ZeroDim_exainple_42 (cf. |G12) ). 

Example 3.6. cpdm5 (cf. fV99] ) 

Table [T] summarizes the results where modSolve(c) denotes the parallelized ver- 
sion of the algorithm applied on c cores. All timings are given in seconds. 



Example 


solve 


modSolve 


modSolve(4) 


modSolve(16) 


13.21 


> Ih 


993.7 


994.0 


701.9 


13.31 


900.9 


84.6 


54.8 


35.0 


13.41 


> 2h 


981.8 


955.2 


846.3 


13.51 


816.0 


643.5 


483.3 


124.3 


13.61 


> Ih 


370.1 


370.0 


243.3 



Table 1 . Total running times in seconds for computing all roots of 
the considered examples via solve, modSolve and its parallelized 
variant modSolve(c) for c = 4, 16. 
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